################################################################################################
##########Daniel Bochsler & Andreas Juon##########
##########Power-sharing and the quality of democracy#########
##########Replication script##########
################################################################################################

#####1. Data import and set-up#####

###data###
setwd("****ENTER WORKING DIRECTORY WHERE REPLICATION FOLDER IS LOCATED****")
psqod <- read.csv("./ps_qod_replication_data.csv", sep=",")
psqod1 <- subset(psqod, country!="Cyprus")
psqod1 <- subset(psqod1, country == "Malta" | country == "Montenegro" | country == "Iceland" | (polity_l1 >= 6 | is.na(polity_l1)) | (polity >= 6 | is.na(polity)))
indicator <- read.csv(file="./indicator_description.csv")

###libraries###
library(tidyr)
library(ggrepel)
library(RColorBrewer)
library(gridExtra)
library(dotwhisker)
library(stargazer)
library(ggplot2)
library(grid)
library(lmtest)
library(multiwayvcov)
library(plyr)
library(dplyr)
library(broom)

###formulas###
get_confint<-function(model, vcovCL){
  t<-qt(.975, model$df.residual)
  ct<-coeftest(model, vcovCL)
  est<-cbind(ct[,1], ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
  colnames(est)<-c("Estimate","LowerCI","UpperCI")
  return(est)
}
cluster.se<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + 
                    theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x +
                 theme(legend.position = "none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl), 
                                            legend,ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend, ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  grid.newpage()
  grid.draw(combined)
  invisible(combined)
}


#####2. Main models#####

###definition of variables###
controls_d1 <- "log(gdppc_l1) + log(pop_l1) + log(country_area_l1) + oilrents_l1 + fractionalization1_l1 + log(nopluralitysum_l1 + 0.0000001) + postconflict_l1 + ongoingc_l1 + region2 + year"#same controls as in IV
principles <- c("INDLIB", "RULEOFLAW", "PUBLIC", "COMPET", "MUTUCONS", "GOVCAP", "TRANSPAR", "PARTICIP", "REPRES")
components <- c("IL_PHIN", "IL_SELFU", "RL_EQL", "RL_QUAL", "PS_FRAS", "PS_FROP", "CO_COMP", "CO_OPEN", "MC_CHECKS", "MC_VERT", "GC_GORE", "GC_CEIM", "TR_NOSEC", "TR_PTPP", "PAR_EQPA", "PAR_EFPA", "REP_SR", "REP_DR")
subcomponents <- c("IL_PHIN1", "IL_PHIN2", "IL_PHIN3", "IL_SELFU1", "IL_SELFU2", "IL_SELFU3", "RL_EQL1", "RL_EQL2", "RL_EQL3", "RL_QUAL1", "RL_QUAL2", "RL_QUAL3", "PS_FRAS1", "PS_FRAS2", "PS_FRAS3", "PS_FROP1", "PS_FROP2", "PS_FROP3", "CO_COMP1", "CO_COMP2", "CO_COMP3", "CO_OPEN1", "CO_OPEN2", "CO_OPEN3", "MC_CHECKS1", "MC_CHECKS2", "MC_CHECKS3", "MC_VERT1", "MC_VERT2", "GC_GORE1", "GC_GORE2", "GC_GORE3", "GC_CEIM1", "GC_CEIM2", "GC_CEIM3", "GC_CEIM4", "TR_NOSEC1", "TR_NOSEC2", "TR_PTPP1", "TR_PTPP2", "TR_PTPP3", "PAR_EQPA1", "PAR_EQPA2", "PAR_EQPA3", "PAR_EFPA1", "PAR_EFPA2", "PAR_EFPA3", "REP_SR1", "REP_SR2", "REP_SR3", "REP_DR1", "REP_DR2", "REP_DR3")
indicators <- c("Consttort", "Convtort", "Politterr", "Torture", "Homicide", "Riot", "Constrel", "Constfreemov", "Freerelig", "Freemove", "Propright", "Secprop", "Constfair", "Pubtrial", "Judindepcor", "Judindepinf", "Impcourts", "Integrlegal", "Profjudg", "Proftenure", "Confjust", "Fairjust", "Confpolice", "Fairpolice", "Constfras", "Constass", "Union", "Memproorg", "Memhuman", "Memenviron", "Constspeech", "Constpress", "Newsimp", "Newspaper", "Balpress", "Neutrnp", "Meandistrict", "Gerryman", "Largpavo", "Votediff", "Herfindex", "Seatdiff", "Adminhurd", "Eff_thresh", "Smallpavo", "Enep", "Ceilings", "Funding", "Balexleg", "Balpowexle", "Seatsgov_2", "Judrev", "Powjudi", "Federalism", "Bicameralism", "Subexp", "Subrev", "Leglen", "Gov_term", "Confgov", "Devbehav", "Govstab", "Cabchange", "Antigovact", "Violantigov", "Mip", "Rip", "Publser", "Govdec", "Bureau", "CenBank_Ind", "Discinco", "Discexp", "Corrup", "CPI", "RestricFoi", "EffFoi", "Legmedia", "Polmedia", "Transp", "Suffrage", "Regprovap", "Repturnined", "Repturngeag", "Repaltined", "Repaltgeag", "Facilitat", "Registr", "Meanpart", "Eff_DD", "Petitions", "Demons", "Seatperinh", "No_districts", "Dirdem", "DD_Quora", "Gallagindex", "Issuecongr", "Polrightwom", "Constraints", "Partreg", "Womrep", "womgov", "Minrep", "Minpower")#volatility

###Model 1###
m1 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m1[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m1)[n] <- i
  n <- n + 1
}
m1_list <- m1
m1 <- ldply(m1, tidy, .id = "model")
m1 <- left_join(m1,indicator,by="model")
m1 <- subset(m1,term == "ps_avg2_l1" & concept_type == "indicator")
m1 <- m1[order(m1$hypothesis_number,m1$concept_type,m1$hypothesis_expectation),]
m1_full <- m1
m1$hypothesis_sig <- ifelse(m1$p.value < 0.1, "*", "")
m1$hypothesis_sig <- ifelse(m1$p.value < 0.05, "**", m1$hypothesis_sig)
m1$hypothesis_sig <- ifelse(m1$p.value < 0.01, "***", m1$hypothesis_sig)
m1$estimate_m1 <- paste(round(m1$estimate,2),m1$hypothesis_sig,sep="")
m1$p.value_m1 <- round(m1$p.value,3)
m1 <- m1[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_m1","p.value_m1")]
write.csv(m1,file="./model_results/m1.csv")


#####3. Robustness checks#####

###model 2 / robustness check 1###
m2 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + pr_elecsys_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m2[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m2)[n] <- i
  n <- n + 1
}
m2_list <- m2
m2 <- ldply(m2, tidy, .id = "model")
m2 <- left_join(m2,indicator,by="model")
m2 <- subset(m2,term == "ps_avg2_l1" & concept_type == "indicator")
m2 <- m2[order(m2$hypothesis_number,m2$concept_type,m2$hypothesis_expectation),]
m2_full <- m2
m2$hypothesis_sig <- ifelse(m2$p.value < 0.1, "*", "")
m2$hypothesis_sig <- ifelse(m2$p.value < 0.05, "**", m2$hypothesis_sig)
m2$hypothesis_sig <- ifelse(m2$p.value < 0.01, "***", m2$hypothesis_sig)
m2$estimate_m2 <- paste(round(m2$estimate,2),m2$hypothesis_sig,sep="")
m2$p.value_m2 <- round(m2$p.value,3)
m2 <- m2[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_m2","p.value_m2")]
write.csv(m2,file="./model_results/m2.csv")

###model 3 / robustness check 2###
m3 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ ps_avg2_nopr_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m3[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m3)[n] <- i
  n <- n + 1
}
m3_list <- m3
m3 <- ldply(m3, tidy, .id = "model")
m3 <- left_join(m3,indicator,by="model")
m3 <- subset(m3,term == "ps_avg2_nopr_l1" & concept_type == "indicator")
m3 <- m3[order(m3$hypothesis_number,m3$concept_type,m3$hypothesis_expectation),]
m3_full <- m3
m3$hypothesis_sig <- ifelse(m3$p.value < 0.1, "*", "")
m3$hypothesis_sig <- ifelse(m3$p.value < 0.05, "**", m3$hypothesis_sig)
m3$hypothesis_sig <- ifelse(m3$p.value < 0.01, "***", m3$hypothesis_sig)
m3$estimate_m3 <- paste(round(m3$estimate,2),m3$hypothesis_sig,sep="")
m3$p.value_m3 <- round(m3$p.value,3)
m3 <- m3[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_m3","p.value_m3")]
write.csv(m3,file="./model_results/m3.csv")

###model 4 / robustness check 3###
m4 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ ps_corp_avg2_l1 + ps_lib_avg2_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m4[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m4)[n] <- i
  n <- n + 1
}
m4_list <- m4
m4 <- ldply(m4, tidy, .id = "model")
m4 <- left_join(m4,indicator,by="model")
m4 <- m4[order(m4$hypothesis_number,m4$concept_type,m4$hypothesis_expectation),]
m4_corp <- subset(m4,term == "ps_corp_avg2_l1" & concept_type == "indicator")
m4_full <- m4
m4_corp$hypothesis_sig <- ifelse(m4_corp$p.value < 0.1, "*", "")
m4_corp$hypothesis_sig <- ifelse(m4_corp$p.value < 0.05, "**", m4_corp$hypothesis_sig)
m4_corp$hypothesis_sig <- ifelse(m4_corp$p.value < 0.01, "***", m4_corp$hypothesis_sig)
m4_corp$estimate_corp_m4 <- paste(round(m4_corp$estimate,2),m4_corp$hypothesis_sig,sep="")
m4_corp$p.value_corp_m4 <- round(m4_corp$p.value,3)
m4_corp <- m4_corp[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_corp_m4","p.value_corp_m4")]
m4_lib <- subset(m4,term == "ps_lib_avg2_l1" & concept_type == "indicator")
m4_lib$hypothesis_sig <- ifelse(m4_lib$p.value < 0.1, "*", "")
m4_lib$hypothesis_sig <- ifelse(m4_lib$p.value < 0.05, "**", m4_lib$hypothesis_sig)
m4_lib$hypothesis_sig <- ifelse(m4_lib$p.value < 0.01, "***", m4_lib$hypothesis_sig)
m4_lib$estimate_lib_m4 <- paste(round(m4_lib$estimate,2),m4_lib$hypothesis_sig,sep="")
m4_lib$p.value_lib_m4 <- round(m4_lib$p.value,3)
m4_lib <- m4_lib[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_lib_m4","p.value_lib_m4")]
m4 <- left_join(m4_corp,m4_lib[c("model","estimate_lib_m4","p.value_lib_m4")], by="model")
write.csv(m4,file="./model_results/m4.csv")

###model 5 / robustness check 4###
psqod1_subset <- subset(psqod1,nogroups1_epr_l1>1)
m5 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, ", data = psqod1_subset)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1_subset$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m5[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m5)[n] <- i
  n <- n + 1
}
m5_list <- m5
m5 <- ldply(m5, tidy, .id = "model")
m5 <- left_join(m5,indicator,by="model")
m5 <- subset(m5,term == "ps_avg2_l1" & concept_type == "indicator")
m5 <- m5[order(m5$hypothesis_number,m5$concept_type,m5$hypothesis_expectation),]
m5_full <- m5
m5$hypothesis_sig <- ifelse(m5$p.value < 0.1, "*", "")
m5$hypothesis_sig <- ifelse(m5$p.value < 0.05, "**", m5$hypothesis_sig)
m5$hypothesis_sig <- ifelse(m5$p.value < 0.01, "***", m5$hypothesis_sig)
m5$estimate_m5 <- paste(round(m5$estimate,2),m5$hypothesis_sig,sep="")
m5$p.value_m5 <- round(m5$p.value,3)
m5 <- m5[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_m5","p.value_m5")]
write.csv(m5,file="./model_results/m5.csv")

###model 6 / robustness check 5###
psqod1_subset <- subset(psqod1,postconflict_l1==1)
m6 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, ", data = psqod1_subset)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1_subset$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m6[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m6)[n] <- i
  n <- n + 1
}
m6_list <- m6
m6 <- ldply(m6, tidy, .id = "model")
m6 <- left_join(m6,indicator,by="model")
m6 <- subset(m6,term == "ps_avg2_l1" & concept_type == "indicator")
m6 <- m6[order(m6$hypothesis_number,m6$concept_type,m6$hypothesis_expectation),]
m6_full <- m6
m6$hypothesis_sig <- ifelse(m6$p.value < 0.1, "*", "")
m6$hypothesis_sig <- ifelse(m6$p.value < 0.05, "**", m6$hypothesis_sig)
m6$hypothesis_sig <- ifelse(m6$p.value < 0.01, "***", m6$hypothesis_sig)
m6$estimate_m6 <- paste(round(m6$estimate,2),m6$hypothesis_sig,sep="")
m6$p.value_m6 <- round(m6$p.value,3)
m6 <- m6[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_m6","p.value_m6")]
write.csv(m6,file="./model_results/m6.csv")

###model 7 / robustness check 6###
m7 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, "+ regime_durability", ", data = psqod1)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m7[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m7)[n] <- i
  n <- n + 1
}
m7_list <- m7
m7 <- ldply(m7, tidy, .id = "model")
m7 <- left_join(m7,indicator,by="model")
m7 <- subset(m7,term == "ps_avg2_l1" & concept_type == "indicator")
m7 <- m7[order(m7$hypothesis_number,m7$concept_type,m7$hypothesis_expectation),]
m7_full <- m7
m7$hypothesis_sig <- ifelse(m7$p.value < 0.1, "*", "")
m7$hypothesis_sig <- ifelse(m7$p.value < 0.05, "**", m7$hypothesis_sig)
m7$hypothesis_sig <- ifelse(m7$p.value < 0.01, "***", m7$hypothesis_sig)
m7$estimate_m7 <- paste(round(m7$estimate,2),m7$hypothesis_sig,sep="")
m7$p.value_m7 <- round(m7$p.value,3)
m7 <- m7[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_m7","p.value_m7")]
write.csv(m7,file="./model_results/m7.csv")

###model 8 / robustness check 7###
m8 <- list()
n <- 1
for (i in (c(principles, components, subcomponents, indicators))) {
  step1 <- paste(i,"_unclustered <- lm(", i, " ~ inclusive_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste(i," <- cluster.se(",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  step3 <- paste("rm(",i,"_unclustered)",sep="")
  step4 <- paste("m8[[",n,"]] <- ",i,sep="")
  step5 <- paste("rm(",i,")",sep="")
  eval(parse(text=c(step1, step2, step3, step4, step5)))
  names(m8)[n] <- i
  n <- n + 1
}
m8_list <- m8
m8 <- ldply(m8, tidy, .id = "model")
m8 <- left_join(m8,indicator,by="model")
m8 <- subset(m8,term == "inclusive_l1" & concept_type == "indicator")
m8 <- m8[order(m8$hypothesis_number,m8$concept_type,m8$hypothesis_expectation),]
m8_full <- m8
m8$hypothesis_sig <- ifelse(m8$p.value < 0.1, "*", "")
m8$hypothesis_sig <- ifelse(m8$p.value < 0.05, "**", m8$hypothesis_sig)
m8$hypothesis_sig <- ifelse(m8$p.value < 0.01, "***", m8$hypothesis_sig)
m8$estimate_m8 <- paste(round(m8$estimate,2),m8$hypothesis_sig,sep="")
m8$p.value_m8 <- round(m8$p.value,3)
m8 <- m8[c("model","short_description","layer","hypothesis_number","hypothesis_expectation","estimate_m8","p.value_m8")]
write.csv(m8,file="./model_results/m8.csv")


#####4. Figures and tables in main article + their expanded versions in appendix#####

###table 1###
table1 <- psqod1[c("country","year","ps_avg2")]
table1$psh1_strength4_period <- round(as.numeric(table1$ps_avg2), digits = 1)
table1_lag <- table1
table1_lag$year <- table1_lag$year + 1
table1_lag$psh1_strength4_period_lag <- table1_lag$psh1_strength4_period
table1 <- left_join(table1,table1_lag[c("country","year","psh1_strength4_period_lag")],by=c("country","year"))
table1$psh1_strength4_period_lag <- ifelse(is.na(table1$psh1_strength4_period_lag),-999,table1$psh1_strength4_period_lag)
table1$psh1_strength4_period <- ifelse(is.na(table1$psh1_strength4_period),-999,table1$psh1_strength4_period)
table1$pschange <- ifelse(table1$psh1_strength4_period != table1$psh1_strength4_period_lag,1,0)
table1$pschange <- ifelse(is.na(table1$pschange),1,table1$pschange)
table1 <- table1 %>% group_by(country) %>% arrange(country,year) %>% mutate(period = cumsum(pschange))
table1$psh1_strength4_period <- ifelse(table1$psh1_strength4_period == -999, NA, table1$psh1_strength4_period)
table1 <- table1 %>% group_by(country,period) %>% arrange(country,period) %>% mutate(from = min(year))
table1 <- table1 %>% group_by(country,period) %>% arrange(country,period) %>% mutate(to = max(year))
table1 <- table1 %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_psh = min(ps_avg2))
table1$min_psh <- round(as.numeric(table1$min_psh),digits=2)
table1 <- table1 %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_psh = max(ps_avg2))
table1$max_psh <- round(as.numeric(table1$max_psh),digits=2)
table1$time_period <- ifelse(table1$from != table1$to, paste(table1$from,"-",table1$to,sep=""),table1$from)
table1$horizontal_PS <- ifelse(table1$min_psh != table1$max_psh, paste(table1$min_psh,"-",table1$max_psh,sep=""),table1$max_psh)
table1 <- table1[order(-table1$max_psh),]
table1 <- unique(table1[c("country","time_period","horizontal_PS")])
table1 <- subset(table1,!is.na(horizontal_PS))
write.csv(table1, file="./main_figures/table1_full.csv")

###figure 2###
figure2_data <- psqod[c("country","year","Minrep","Womrep","Constrel","Rip","Govstab","Antigovact","Govdec","Subexp","Judindepinf","CPI")]
figure2_data <- subset(figure2_data, country!= "Cyprus" & !(country == "Kosovo" & year < 2008))
##calculate avg for other countries
figure2_data$country <- as.factor(ifelse(figure2_data$country != "Bosnia-Herzegovina" & figure2_data$country != "Czechoslovakia" & figure2_data$country != "Belgium" & figure2_data$country != "South Africa" & figure2_data$country != "Kosovo" & figure2_data$country != "Switzerland" & figure2_data$country != "Macedonia", "non power-sharing", paste(figure2_data$country)))
for (i in (c("Minrep","Womrep","Constrel","Rip","Govstab","Antigovact","Govdec","Subexp","Judindepinf","CPI"))) {
  step1 <- paste("figure2_data <- figure2_data %>% group_by(country,year)%>% mutate(",i," = mean(",i,",na.rm=T))",sep="")
  eval(parse(text=c(step1)))
}
figure2_data <- unique(figure2_data)
colours <- brewer.pal(8, "Set1")
colours[6] <- "#000000"
##Minrep
figure2_Minrep <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Minrep, group = country, colour=country, linetype =country, size = country)) + ylab("Minrep") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Womrep
figure2_Womrep <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Womrep, group = country, colour=country, linetype =country, size = country)) + ylab("Womrep") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Constrel
figure2_Constrel <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Constrel, group = country, colour=country, linetype =country, size = country)) + ylab("Constrel") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) +scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Rip
figure2_Rip <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Rip, group = country, colour=country, linetype =country, size = country)) + ylab("Rip") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Govstab
figure2_Govstab <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Govstab, group = country, colour=country, linetype =country, size = country)) + ylab("Govstab") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Antigovact
figure2_Antigovact <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Antigovact, group = country, colour=country, linetype =country, size = country)) + ylab("Antigovact") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Govdec
figure2_Govdec <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Govdec, group = country, colour=country, linetype =country, size = country)) + ylab("Govdec") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Subexp
figure2_Subexp <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Subexp, group = country, colour=country, linetype =country, size = country)) + ylab("Subexp") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##Judindepinf
figure2_Judindepinf <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=Judindepinf, group = country, colour=country, linetype =country, size = country)) + ylab("Judindepinf") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##CPI
figure2_CPI <- ggplot(data = figure2_data) + geom_line(aes(x=year, y=CPI, group = country, colour=country, linetype =country, size = country)) + ylab("CPI") + xlab("Year") + scale_color_manual(name = "Country", values = colours) + scale_linetype_manual(name = "Country", values = c(1,2,3,4,5,6,7,8)) + scale_size_manual(name = "Country", values = c(0.5,0.5,0.5,0.5,0.5,0.33,0.5,0.5)) + theme_bw() + theme(legend.position = "bottom", text=element_text(family="Times")) + xlim(1990,2016)+ coord_cartesian(ylim=c(-100, 100))
##overall
figure2 <- grid_arrange_shared_legend(figure2_Minrep + ggtitle("a) Ethnic groups in parliament"), figure2_Womrep + ggtitle("b) Women in parliament"), figure2_Constrel + ggtitle("c) Formal religious freedom"), figure2_Rip + ggtitle("d) Religious interference in politics"), figure2_Govstab + ggtitle("e) Government stability"), figure2_Antigovact + ggtitle("f) Demonstrations and strikes"), figure2_CPI + ggtitle("g) Perception of corruption"), figure2_Subexp + ggtitle("h) Subnational expenditure"), figure2_Judindepinf + ggtitle("i) Judicial independence"), figure2_Govdec + ggtitle("j) Implementation of gov. decisions"), ncol=2,nrow=5)
ggsave(figure2, file="./main_figures/figure2.png", width = 17, height = 25.7, units ="cm", dpi = 300)

###table 2 / table a.8/a.9 [full - main model and robustness checks]###
table2 <- left_join(m1, m2[c("model","estimate_m2","p.value_m2")], by="model")
table2 <- left_join(table2, m3[c("model","estimate_m3","p.value_m3")], by="model")
table2 <- left_join(table2, m4[c("model","estimate_corp_m4","p.value_corp_m4","estimate_lib_m4","p.value_lib_m4")], by="model")
table2 <- left_join(table2, m5[c("model","estimate_m5","p.value_m5")], by="model")
table2 <- left_join(table2, m6[c("model","estimate_m6","p.value_m6")], by="model")
table2 <- left_join(table2, m7[c("model","estimate_m7","p.value_m7")], by="model")
table2 <- left_join(table2, m8[c("model","estimate_m8","p.value_m8")], by="model")
table2[c("p.value_m1","p.value_m2","p.value_m3","p.value_corp_m4","p.value_lib_m4","p.value_m5","p.value_m6","p.value_m7","p.value_m8")] <- NULL
write.csv(table2,file="./main_figures/table2_full.csv")

###figure 3###
##set up dfs
m4_full_corp <- subset(m4_full,term == "ps_corp_avg2_l1" & concept_type == "indicator")
m4_full_lib <- subset(m4_full,term == "ps_lib_avg2_l1" & concept_type == "indicator")
m1_full$model_description <- "M1 (all countries)"
m1_full$model_ps <- "Overall"
m2_full$model_description <- "M2 (all countries, PR control)"
m2_full$model_ps <- "Overall"
m3_full$model_description <- "M3 (all countries, no PR)"
m3_full$model_ps <- "Overall"
m4_full_corp$model_description <- "M4 (all countries, corporate)"
m4_full_corp$model_ps <- "Corporate"
m4_full_lib$model_description <- "M4 (all countries, liberal)"
m4_full_lib$model_ps <- "Liberal"
m5_full$model_description <- "M5 (multi-ethnic countries)"
m5_full$model_ps <- "Overall"
m6_full$model_description <- "M6 (post-conflict countries)"
m6_full$model_ps <- "Overall"
m7_full$model_description <- "M7 (all countries)"
m7_full$model_ps <- "Overall"
m8_full$model_description <- "M8 (all countries)"
m8_full$model_ps <- "Inclusive"
for (i in (c("m1_full","m2_full","m3_full","m4_full_corp","m4_full_lib","m5_full","m6_full","m7_full","m8_full"))) {
  step1 <- paste("", i, "$conf.low <- ", i, "$estimate - 1.645 *", i, "$std.error",sep="")
  step2 <- paste("", i, "$conf.high <- ", i, "$estimate + 1.645 *", i, "$std.error",sep="")
  step3 <- paste("", i, "$model2 <- ", i, "$term",sep="")
  eval(parse(text=c(step1, step2, step3)))
}
figure3_df <- rbind.fill(m1_full, m2_full, m3_full, m4_full_corp, m4_full_lib, m5_full, m6_full, m7_full, m8_full)
figure3_df <- subset(figure3_df, model == "Minrep" | model == "Womrep" | model == "Constrel" | model == "Rip" | model == "Govstab" | model == "Antigovact" | model == "Govdec" | model == "Subexp" | model == "Judindepinf" | model == "CPI")
figure3_df$term <- as.factor(figure3_df$model)
figure3_df$term = factor(figure3_df$term,levels(figure3_df$term)[c(7,10,2,8,5,1,4,9,6,3)])
figure3_df$model <- as.factor(figure3_df$model_description)
figure3_df$group_ps <- paste(figure3_df$model_description, figure3_df$model_ps, sep="_")
figure3_df$group_ps <- as.factor(figure3_df$group_ps)
figure3_df$group_ps = factor(figure3_df$group_ps,levels(figure3_df$group_ps)[c(7,6,5,4,3,2,1)])
figure3_df$model_ps <- as.factor(figure3_df$model_ps)
figure3_df$model_ps = factor(figure3_df$model_ps,levels(figure3_df$model_ps)[c(3,1,2)])
figure3_df$model_description <- as.factor(figure3_df$model_description)
figure3_df$model_description = factor(figure3_df$model_description,levels(figure3_df$model_description)[c(9,8,7,6,5,4,3,2,1)])
figure3_df <- figure3_df%>% arrange(desc(model_description))
figure3_df$xb <- 0
figure3_df$yb<-0
##figure 3 full (appendix)
figure3_full <- dwplot(figure3_df, dodge_size = 0.8,
                       vline = geom_vline(xintercept = 0, colour = "black", linetype = 2),
                       dot_args = list(aes(group = model_description, colour = factor(model_description), shape = factor(model_description))),
                       whisker_args = list(aes(group = model_description, colour = factor(model_description), linetype = factor(model_description)))) +
  theme_bw() + xlab("Coefficient Estimate") + ylab("") +
  geom_line(data = figure3_df, aes(x= xb, y= yb, linetype = factor(model_description))) +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  scale_color_manual(name = "Model", values = c("black","black","grey70","grey70","black","black","grey45","grey45","black"), breaks = levels(as.factor(figure3_df$model_description))) +
  scale_shape_manual(name = "Model", values = c(6,8,17,15,16,16,17,15,16), breaks = levels(as.factor(figure3_df$model_description))) +
  scale_linetype_manual(name = "Model", values = c(1,1,1,1,3,2,1,1,1), breaks = levels(as.factor(figure3_df$model_description)), na.translate=FALSE) + 
  theme(legend.position = "bottom", text=element_text(family="Times")) +guides(colour=guide_legend(reverse = TRUE, title.position="top", nrow=3,byrow=TRUE),shape=guide_legend(reverse = TRUE, title.position="top", nrow=3,byrow=TRUE),linetype=guide_legend(reverse = TRUE, title.position="top", nrow=2,byrow=TRUE))
ggsave(file="./appendix_figures/figure3_full.png", figure3_full, width = 20, height = 20, units ="cm", dpi = 300)
##figure 3 (main article)
figure3_df_small <- subset(figure3_df, model_description == "M1 (all countries)" | model_description == "M4 (all countries, liberal)" | model_description == "M4 (all countries, corporate)")
figure3_small <- dwplot(figure3_df_small, dodge_size = 0.8,
                        vline = geom_vline(xintercept = 0, colour = "black", linetype = 2),
                        dot_args = list(aes(group = model_description, colour = factor(model_description), shape = factor(model_description))),
                        whisker_args = list(aes(group = model_description, colour = factor(model_description), linetype = factor(model_description)))) +
  theme_bw() + xlab("Coefficient Estimate") + ylab("") +
  geom_line(data = figure3_df_small, aes(x= xb, y= yb, linetype = factor(model_description))) +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  scale_color_manual(name = "Model", values = c("black","black","black")) +
  scale_shape_manual(name = "Model", values = c(16,16,16)) +
  scale_linetype_manual(name = "Model", values = c(3,2,1), na.translate=FALSE) + 
  theme(legend.position = "bottom", text=element_text(family="Times")) +guides(colour=guide_legend(reverse = TRUE, title.position="top", nrow=2,byrow=TRUE),shape=guide_legend(reverse = TRUE, title.position="top", nrow=2,byrow=TRUE),linetype=guide_legend(reverse = TRUE, title.position="top", nrow=2,byrow=TRUE))
ggsave(file="./main_figures/figure3.png", figure3_small, width = 20, height = 12, units ="cm", dpi = 300)

###figure 4###
figure4_data <- psqod[c("country","year","Minrep","Womrep","Constrel","Rip","Govstab","Antigovact","Govdec","Subexp","Judindepinf","CPI")]
figure4_data <- gather(figure4_data, variable, outcome, Minrep:CPI, factor_key=TRUE)
figure4 <- ggplot(data = subset(figure4_data, country == "Macedonia" & !is.na(outcome)), aes(x=year, y=outcome)) + geom_line(aes(x=year, y=outcome, group = variable, colour=variable, linetype =variable, size = variable)) +
  ylab("Value of indicator") + xlab("Year") +
  scale_color_manual(name = "Indicator", values = c("black","black","grey45","black","grey70","grey45","grey45","grey45")) +
  scale_linetype_manual(name = "Indicator", values = c(1,2,1,3,1,2,3,4)) +
  scale_size_manual(name = "Indicator", values = c(0.75,0.5,0.75,0.5,0.5,0.75,0.75,0.5), na.translate=FALSE, guide = guide_legend(reverse = TRUE)) + 
  theme_bw() + theme(legend.position = "none", text=element_text(family="Times")) +
  xlim(1990,2022)+
  geom_label_repel(data = subset(figure4_data, country == "Macedonia" & !is.na(outcome) & year==2016), aes(family="Times", label = variable),
                   box.padding   = 0.25, 
                   point.padding = 0.25,
                   size = 3, nudge_y = 20, nudge_x = 5, segment.alpha = 0.5,
                   segment.color = 'grey50')
ggsave(figure4, file="./main_figures/figure4.png", width = 17, height = 8, units ="cm", dpi = 300)


#####5. Other appendix figures and tables#####

###table a.7###
tablea7_data <- psqod1[c("ps_avg2_l1","ps_corp_avg2","ps_lib_avg2","gdppc_l1", "pop_l1", "country_area_l1", "oilrents_l1", "fractionalization1_l1", "nopluralitysum_l1", "postconflict_l1", "ongoingc_l1", "region2", "year",indicators)]
tablea7_data <- data.frame(tablea7_data)
stargazer(tablea7_data, type ="text",title = "Table A7. Descriptive statistics")

###tables a.10 - a.16###
for (i in (c(indicators))) {
  step1 <- paste("m1_",i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste("m1_",i,"_cse <- cluster.se(m1_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m1_Minrep_unclustered, m1_Minpower_unclustered, m1_Partreg_unclustered, m1_Repturnined_unclustered, m1_Repturngeag_unclustered, m1_Repaltined_unclustered, m1_Repaltgeag_unclustered, m1_Issuecongr_unclustered, m1_Polrightwom_unclustered, m1_Womrep_unclustered, m1_womgov_unclustered, se = c(data.frame(m1_Minrep_cse[, 2]), data.frame(m1_Minpower_cse[, 2]), data.frame(m1_Partreg_cse[, 2]), data.frame(m1_Repturnined_cse[, 2]), data.frame(m1_Repturngeag_cse[, 2]), data.frame(m1_Repaltined_cse[, 2]), data.frame(m1_Repaltgeag_cse[, 2]), data.frame(m1_Issuecongr_cse[, 2]), data.frame(m1_Polrightwom_cse[, 2]), data.frame(m1_Womrep_cse[, 2]), data.frame(m1_womgov_cse[, 2])), type ="text",title = "Table A10. Full results: M1 (indicators belonging to hypotheses 1-2)", covariate.labels = c("Power-sharing","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))
stargazer(m1_Consttort_unclustered, m1_Convtort_unclustered, m1_Torture_unclustered, m1_Constrel_unclustered, m1_Freerelig_unclustered, m1_Constfras_unclustered, m1_Constass_unclustered, m1_Constspeech_unclustered, m1_Constpress_unclustered, m1_Constfreemov_unclustered, m1_Freemove_unclustered, m1_Rip_unclustered, se = c(data.frame(m1_Consttort_cse[, 2]), data.frame(m1_Convtort_cse[, 2]), data.frame(m1_Torture_cse[, 2]), data.frame(m1_Constrel_cse[, 2]), data.frame(m1_Freerelig_cse[, 2]), data.frame(m1_Constfras_cse[, 2]), data.frame(m1_Constass_cse[, 2]), data.frame(m1_Constspeech_cse[, 2]), data.frame(m1_Constpress_cse[, 2]), data.frame(m1_Constfreemov_cse[, 2]), data.frame(m1_Freemove_cse[, 2]), data.frame(m1_Rip_cse[, 2])), type ="text",title = "Table A11. Full results: M1 (indicators belonging to hypotheses 3-4)", covariate.labels = c("Power-sharing","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))
stargazer(m1_Balpowexle_unclustered, m1_Cabchange_unclustered, m1_Largpavo_unclustered, m1_Votediff_unclustered, m1_Herfindex_unclustered, m1_Seatdiff_unclustered, m1_Smallpavo_unclustered, m1_Enep_unclustered, m1_Seatsgov_2_unclustered, m1_Govstab_unclustered, m1_Gov_term_unclustered, m1_Leglen_unclustered, m1_Gerryman_unclustered, se = c(data.frame(m1_Balpowexle_cse[, 2]), data.frame(m1_Cabchange_cse[, 2]), data.frame(m1_Largpavo_cse[, 2]), data.frame(m1_Votediff_cse[, 2]), data.frame(m1_Herfindex_cse[, 2]), data.frame(m1_Seatdiff_cse[, 2]), data.frame(m1_Smallpavo_cse[, 2]), data.frame(m1_Enep_cse[, 2]), data.frame(m1_Seatsgov_2_cse[, 2]), data.frame(m1_Govstab_cse[, 2]), data.frame(m1_Gov_term_cse[, 2]), data.frame(m1_Leglen_cse[, 2]), data.frame(m1_Gerryman_cse[, 2])), type ="text",title = "Table A12: Full results: M1 (indicators belonging to hypothesis 5)", covariate.labels = c("Power-sharing","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))
stargazer(m1_Union_unclustered, m1_Memproorg_unclustered, m1_Memhuman_unclustered, m1_Memenviron_unclustered, m1_Antigovact_unclustered, m1_Regprovap_unclustered, m1_Meanpart_unclustered, m1_Petitions_unclustered, m1_Demons_unclustered, se = c(data.frame(m1_Union_cse[, 2]), data.frame(m1_Memproorg_cse[, 2]), data.frame(m1_Memhuman_cse[, 2]), data.frame(m1_Memenviron_cse[, 2]), data.frame(m1_Antigovact_cse[, 2]), data.frame(m1_Regprovap_cse[, 2]), data.frame(m1_Meanpart_cse[, 2]), data.frame(m1_Petitions_cse[, 2]), data.frame(m1_Demons_cse[, 2])), type ="text",title = "Table A13. Full results: M1 (indicators belonging to hypothesis 6)", covariate.labels = c("Power-sharing","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))
stargazer(m1_Ceilings_unclustered, m1_Funding_unclustered, m1_Publser_unclustered, m1_Bureau_unclustered, m1_Discinco_unclustered, m1_Discexp_unclustered, m1_Corrup_unclustered, m1_CPI_unclustered, m1_RestricFoi_unclustered, m1_EffFoi_unclustered, m1_Transp_unclustered, se = c(data.frame(m1_Ceilings_cse[, 2]), data.frame(m1_Funding_cse[, 2]), data.frame(m1_Publser_cse[, 2]), data.frame(m1_Bureau_cse[, 2]), data.frame(m1_Discinco_cse[, 2]), data.frame(m1_Discexp_cse[, 2]), data.frame(m1_Corrup_cse[, 2]), data.frame(m1_CPI_cse[, 2]), data.frame(m1_RestricFoi_cse[, 2]), data.frame(m1_EffFoi_cse[, 2]), data.frame(m1_Transp_cse[, 2])), type ="text",title = "Table A14. Full results: M1 (indicators belonging to hypothesis 7)", covariate.labels = c("Power-sharing","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))
stargazer(m1_Balexleg_unclustered, m1_Bicameralism_unclustered, m1_Federalism_unclustered, m1_Subexp_unclustered, m1_Subrev_unclustered, m1_Confgov_unclustered, m1_Govdec_unclustered, se = c(data.frame(m1_Balexleg_cse[, 2]), data.frame(m1_Bicameralism_cse[, 2]), data.frame(m1_Federalism_cse[, 2]), data.frame(m1_Subexp_cse[, 2]), data.frame(m1_Subrev_cse[, 2]), data.frame(m1_Confgov_cse[, 2]), data.frame(m1_Govdec_cse[, 2])), type ="text",title = "Table A15. Full results: M1 (indicators belonging to hypotheses 8 and 10)", covariate.labels = c("Power-sharing","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))
stargazer(m1_Judindepcor_unclustered, m1_Judindepinf_unclustered, m1_Impcourts_unclustered, m1_Integrlegal_unclustered, m1_Profjudg_unclustered, m1_Proftenure_unclustered, m1_Confjust_unclustered, m1_Fairjust_unclustered, m1_Judrev_unclustered, m1_Powjudi_unclustered, m1_CenBank_Ind_unclustered, se = c(data.frame(m1_Judindepcor_cse[, 2]), data.frame(m1_Judindepinf_cse[, 2]), data.frame(m1_Impcourts_cse[, 2]), data.frame(m1_Integrlegal_cse[, 2]), data.frame(m1_Profjudg_cse[, 2]), data.frame(m1_Proftenure_cse[, 2]), data.frame(m1_Confjust_cse[, 2]), data.frame(m1_Fairjust_cse[, 2]), data.frame(m1_Judrev_cse[, 2]), data.frame(m1_Powjudi_cse[, 2]), data.frame(m1_CenBank_Ind_cse[, 2])), type ="text",title = "Table A16. Full results: M1 (indicators belonging to hypothesis 9)", covariate.labels = c("Power-sharing","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))

###table a.17###
exemplary_indicators <- c("Minrep","Womrep","Constrel","Rip","Govstab","Antigovact","Govdec","Subexp","Judindepinf","CPI")
for (i in (c(exemplary_indicators))) {
  step1 <- paste("m2_",i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + pr_elecsys_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste("m2_",i,"_cse <- cluster.se(m2_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m2_Minrep_unclustered, m2_Womrep_unclustered, m2_Constrel_unclustered, m2_Rip_unclustered, m2_Govstab_unclustered, m2_Antigovact_unclustered, m2_Govdec_unclustered, m2_Subexp_unclustered, m2_Judindepinf_unclustered, m2_CPI_unclustered, se = c(data.frame(m2_Minrep_cse[, 2]), data.frame(m2_Womrep_cse[, 2]), data.frame(m2_Constrel_cse[, 2]), data.frame(m2_Rip_cse[, 2]), data.frame(m2_Govstab_cse[, 2]), data.frame(m2_Antigovact_cse[, 2]), data.frame(m2_Govdec_cse[, 2]), data.frame(m2_Subexp_cse[, 2]), data.frame(m2_Judindepinf_cse[, 2]), data.frame(m2_CPI_cse[, 2])), type ="text",title = "Table A17. Full results: M2 (robustness check 1, exemplary indicators)")

###table a.18###
for (i in (c(exemplary_indicators))) {
  step1 <- paste("m3_",i,"_unclustered <- lm(", i, " ~ ps_avg2_nopr_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste("m3_",i,"_cse <- cluster.se(m3_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m3_Minrep_unclustered, m3_Womrep_unclustered, m3_Constrel_unclustered, m3_Rip_unclustered, m3_Govstab_unclustered, m3_Antigovact_unclustered, m3_Govdec_unclustered, m3_Subexp_unclustered, m3_Judindepinf_unclustered, m3_CPI_unclustered, se = c(data.frame(m3_Minrep_cse[, 2]), data.frame(m3_Womrep_cse[, 2]), data.frame(m3_Constrel_cse[, 2]), data.frame(m3_Rip_cse[, 2]), data.frame(m3_Govstab_cse[, 2]), data.frame(m3_Antigovact_cse[, 2]), data.frame(m3_Govdec_cse[, 2]), data.frame(m3_Subexp_cse[, 2]), data.frame(m3_Judindepinf_cse[, 2]), data.frame(m3_CPI_cse[, 2])), type ="text",title = "Table A18. Full results: M3 (robustness check 2, exemplary indicators)")

###table a.19###
for (i in (c(exemplary_indicators))) {
  step1 <- paste("m4_",i,"_unclustered <- lm(", i, " ~ ps_corp_avg2_l1 + ps_lib_avg2_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste("m4_",i,"_cse <- cluster.se(m4_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m4_Minrep_unclustered, m4_Womrep_unclustered, m4_Constrel_unclustered, m4_Rip_unclustered, m4_Govstab_unclustered, m4_Antigovact_unclustered, m4_Govdec_unclustered, m4_Subexp_unclustered, m4_Judindepinf_unclustered, m4_CPI_unclustered, se = c(data.frame(m4_Minrep_cse[, 2]), data.frame(m4_Womrep_cse[, 2]), data.frame(m4_Constrel_cse[, 2]), data.frame(m4_Rip_cse[, 2]), data.frame(m4_Govstab_cse[, 2]), data.frame(m4_Antigovact_cse[, 2]), data.frame(m4_Govdec_cse[, 2]), data.frame(m4_Subexp_cse[, 2]), data.frame(m4_Judindepinf_cse[, 2]), data.frame(m4_CPI_cse[, 2])), type ="text",title = "Table A19. Full results: M4 (robustness check 3, exemplary indicators)")

###table a.20###
psqod1_subset <- subset(psqod1,nogroups1_epr_l1>1)
for (i in (c(exemplary_indicators))) {
  step1 <- paste("m5_",i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, ", data = psqod1_subset)",sep="")
  step2 <- paste("m5_",i,"_cse <- cluster.se(m5_",i,"_unclustered,as.integer(psqod1_subset$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m5_Minrep_unclustered, m5_Womrep_unclustered, m5_Constrel_unclustered, m5_Rip_unclustered, m5_Govstab_unclustered, m5_Antigovact_unclustered, m5_Govdec_unclustered, m5_Subexp_unclustered, m5_Judindepinf_unclustered, m5_CPI_unclustered, se = c(data.frame(m5_Minrep_cse[, 2]), data.frame(m5_Womrep_cse[, 2]), data.frame(m5_Constrel_cse[, 2]), data.frame(m5_Rip_cse[, 2]), data.frame(m5_Govstab_cse[, 2]), data.frame(m5_Antigovact_cse[, 2]), data.frame(m5_Govdec_cse[, 2]), data.frame(m5_Subexp_cse[, 2]), data.frame(m5_Judindepinf_cse[, 2]), data.frame(m5_CPI_cse[, 2])), type ="text",title = "Table A20. Full results: M5 (robustness check 4, exemplary indicators)")

###table a.21###
psqod1_subset <- subset(psqod1,postconflict_l1==1)
for (i in (c(exemplary_indicators))) {
  step1 <- paste("m6_",i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, ", data = psqod1_subset)",sep="")
  step2 <- paste("m6_",i,"_cse <- cluster.se(m6_",i,"_unclustered,as.integer(psqod1_subset$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m6_Minrep_unclustered, m6_Womrep_unclustered, m6_Constrel_unclustered, m6_Rip_unclustered, m6_Govstab_unclustered, m6_Antigovact_unclustered, m6_Govdec_unclustered, m6_Subexp_unclustered, m6_Judindepinf_unclustered, m6_CPI_unclustered, se = c(data.frame(m6_Minrep_cse[, 2]), data.frame(m6_Womrep_cse[, 2]), data.frame(m6_Constrel_cse[, 2]), data.frame(m6_Rip_cse[, 2]), data.frame(m6_Govstab_cse[, 2]), data.frame(m6_Antigovact_cse[, 2]), data.frame(m6_Govdec_cse[, 2]), data.frame(m6_Subexp_cse[, 2]), data.frame(m6_Judindepinf_cse[, 2]), data.frame(m6_CPI_cse[, 2])), type ="text",title = "Table A21. Full results: M6 (robustness check 5, exemplary indicators)")

###table a.22###
for (i in (c(exemplary_indicators))) {
  step1 <- paste("m7_",i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, "+ regime_durability", ", data = psqod1)",sep="")
  step2 <- paste("m7_",i,"_cse <- cluster.se(m7_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m7_Minrep_unclustered, m7_Womrep_unclustered, m7_Constrel_unclustered, m7_Rip_unclustered, m7_Govstab_unclustered, m7_Antigovact_unclustered, m7_Govdec_unclustered, m7_Subexp_unclustered, m7_Judindepinf_unclustered, m7_CPI_unclustered, se = c(data.frame(m7_Minrep_cse[, 2]), data.frame(m7_Womrep_cse[, 2]), data.frame(m7_Constrel_cse[, 2]), data.frame(m7_Rip_cse[, 2]), data.frame(m7_Govstab_cse[, 2]), data.frame(m7_Antigovact_cse[, 2]), data.frame(m7_Govdec_cse[, 2]), data.frame(m7_Subexp_cse[, 2]), data.frame(m7_Judindepinf_cse[, 2]), data.frame(m7_CPI_cse[, 2])), type ="text",title = "Table A22. Full results: m7 (robustness check 5, exemplary indicators)")

###table a.23###
for (i in (c(exemplary_indicators))) {
  step1 <- paste("m8_",i,"_unclustered <- lm(", i, " ~ inclusive_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste("m8_",i,"_cse <- cluster.se(m8_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m8_Minrep_unclustered, m8_Womrep_unclustered, m8_Constrel_unclustered, m8_Rip_unclustered, m8_Govstab_unclustered, m8_Antigovact_unclustered, m8_Govdec_unclustered, m8_Subexp_unclustered, m8_Judindepinf_unclustered, m8_CPI_unclustered, se = c(data.frame(m8_Minrep_cse[, 2]), data.frame(m8_Womrep_cse[, 2]), data.frame(m8_Constrel_cse[, 2]), data.frame(m8_Rip_cse[, 2]), data.frame(m8_Govstab_cse[, 2]), data.frame(m8_Antigovact_cse[, 2]), data.frame(m8_Govdec_cse[, 2]), data.frame(m8_Subexp_cse[, 2]), data.frame(m8_Judindepinf_cse[, 2]), data.frame(m8_CPI_cse[, 2])), type ="text",title = "Table A23. Full results: m8 (robustness check 5, exemplary indicators)")

###table a.24###
for (i in (c("Minrep","Minpower","Partreg","Repturnined","Repturngeag","Repaltined","Repaltgeag","Issuecongr","Polrightwom","Womrep","womgov"))) {
  step1 <- paste("m9_",i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + lh_quota_gender_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste("m9_",i,"_cse <- cluster.se(m9_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m9_Minrep_unclustered, m9_Minpower_unclustered, m9_Partreg_unclustered, m9_Repturnined_unclustered, m9_Repturngeag_unclustered, m9_Repaltined_unclustered, m9_Repaltgeag_unclustered, m9_Issuecongr_unclustered, m9_Polrightwom_unclustered, m9_Womrep_unclustered, m9_womgov_unclustered, se = c(data.frame(m9_Minrep_cse[, 2]), data.frame(m9_Minpower_cse[, 2]), data.frame(m9_Partreg_cse[, 2]), data.frame(m9_Repturnined_cse[, 2]), data.frame(m9_Repturngeag_cse[, 2]), data.frame(m9_Repaltined_cse[, 2]), data.frame(m9_Repaltgeag_cse[, 2]), data.frame(m9_Issuecongr_cse[, 2]), data.frame(m9_Polrightwom_cse[, 2]), data.frame(m9_Womrep_cse[, 2]), data.frame(m9_womgov_cse[, 2])), type ="text",title = "Table A24. Full results: M9 (indicators belonging to hypotheses 1-2)", covariate.labels = c("Power-sharing","Gender quota","GDP pc. (log)","Population (log)", "Area (log)","Fuel exports","Ethnic fractionalization","Minority population (log)", "Post-conflict","Ongoing conflict","Latin America", "North America","Rest of the World","Western Europe", "Year"))

###table a.25###
tablea25_data <- subset(psqod1, lh_quota_gender == 1)
tablea25_data <- tablea25_data %>% group_by(country)  %>% mutate(from = min(year))
tablea25_data <- tablea25_data %>% group_by(country)  %>% mutate(to = max(year))
tablea25_data <- unique(tablea25_data[c("country","from","to")])
write.csv(tablea25_data, file="./appendix_figures/tablea25.csv")

###table a.27###
for (i in (c("v2x_polyarchy","v2x_libdem", "v2x_partipdem", "v2x_delibdem","v2x_egaldem"))) {
  step1 <- paste("m10_",i,"_unclustered <- lm(", i, " ~ ps_avg2_l1 + ", controls_d1, ", data = psqod1)",sep="")
  step2 <- paste("m10_",i,"_cse <- cluster.se(m10_",i,"_unclustered,as.integer(psqod1$gwid))",sep="")
  eval(parse(text=c(step1, step2)))
}
stargazer(m10_v2x_polyarchy_unclustered, m10_v2x_libdem_unclustered, m10_v2x_partipdem_unclustered, m10_v2x_delibdem_unclustered, m10_v2x_egaldem_unclustered, se = c(data.frame(m10_v2x_polyarchy_cse[, 2]), data.frame(m10_v2x_libdem_cse[, 2]), data.frame(m10_v2x_partipdem_cse[, 2]), data.frame(m10_v2x_delibdem_cse[, 2]), data.frame(m10_v2x_egaldem_cse[, 2])), type ="text",title = "Table A27. Full results: V-Dem (alternative dependent variables)")

